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The BCS and/or HFB theories are extended by treating the effect of four quasi-particle states 
perturbatively. The approach is tested on the pairing hamiltonian, showing that it combines the 
advantage of standard perturbation theory valid at low pairing strength and of non-perturbative 
approaches breaking particle number valid at higher pairing strength. Including the restoration of 
particle number, further improves the description of pairing correlation. In the presented test, the 
agreement between the exact solution and the combined perturbative + projection is almost perfect. 
The proposed method scales friendly when the number of particles increases and provides a simple 
alternative to other more complicated approaches. 
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Theories breaking explicitly the particle number, 
like Bardeen-Cooper-Schrieffer (BCS) or Hartree-Fock- 
Bogoliubov (HFB) are the simplest efficient ways to de- 
scribe pairing correlations in physical systems. How- 
ever, these theories suffer from some limitations espe- 
cially when the number of particles becomes rather small 
as it is the case in mesoscopic systems like in nuclear 
physics [l| or condensed matter The main difficulty 
is the sharp transition from normal to superfluid phases 
as the pairing strength increases. The second difficulty is 
a systematic underestimation of pairing correlations that 
increases when the particle number decreases. 

In the last decades, several approaches have been used 
to overcome these difficulties. When the number of parti- 
cle is small enough, the exact solution of the pairing prob- 
lem is accessible by direct diagonalization of the Hamil- 
tonian 0-13 and/or using the secular equation originally 
proposed by Richardson Accurate description of 

the ground state energy of superfluid systems with larger 
particle number can be obtained using either Quantum 
Monte-Carlo technique (see for instance 0, [l(|) or ex- 
tending quasi-particle theories using Variation After Pro- 
jection approaches [HI [l2| (for recent applications see 
[13l - [l5l ]). All of these methods are however rather in- 
volved and demanding in terms of computational power. 
The aim of the present work is to show that a perturba- 
tive approach can eventually provide a simpler and yet 
accurate alternative to treat the pairing problem. 

In the following, standard perturbation theory is first 
recalled. It is then explained how to take advantage 
of both perturbative approaches and BCS/HFB theory. 
The importance of particle number restoration is finally 
highlighted. 



Our starting point here is similar to the one used in 
refs. [H, [l(|. We assume a two-body pairing hamiltonian 
given by: 



H = 



Q 

E 

i=l 



Ei(a\a % + atcq) + ^ v ij a i a q a ] a j (!) 
H + V. 



Here, (i, i) denotes time-reversed pair. Such an hamil- 
tonian can be for instance constructed starting from a 
realistic self-consistent mean-field calculation providing 
single- particle energies and two-body matrix elements 0] . 

Here fl denotes the maximal number of accessible 
levels for the pairs of particles. In the following, we 
will consider the specific case of equidistant levels with 
£j = (i — l)Ae with Ae = 1 McV and a number of parti- 
cles N — CI. 

As already noted some times ago [13] and recently re- 
discussed 0, EH , the two-body part V can be treated 
perturbatively to provide an accurate description of the 
pairing problem in the weak pairing interaction regime. 
Starting from the ground state |$o) of Hq, that corre- 
sponds to the Slater determinant obtained by occupying 
the lowest single-particle states while other excited states 
|$„) of H can be obtained by considering particle-hole 
(p-h), 2p-2h, ... excitations built on top of the |$o)- Us- 
ing second-order perturbation theory, the ground state 
energy of H reads: 



Eq 



E + E t 



(2) 



(2) 



(with iVp air = N/2) is the 



where E = 2^- 

ground state energy of Hq while E^ 1 is the standard 
second-order correction, that can be found in textbook: 



v pair 
(2) 
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In the specific case considered here, i.e. Eq. ([I), only 
2p-2h states contribute to the second order correction. 
More precisely, the states n could only differ from the 
ground state by one occupied pair above the Fermi sea 
(denoted by j) and one unoccupied pair below (denoted 
by i): = ajatajaj|$o) and are eigenstates of Hq 

with energy = Eq + 2(Ej — £j). This leads to 



E, 



(2) 



\ E E 

i=l,W pair 3=JV paJt +l,JV n 



(4) 



with A?p a i r = N/2 and N m! rx = 0/2. This expression 
has been obtained in ref. [18| using a different approach 
based on the Richardson- Gaudin equation. It is known 
fl6l [I?} that standard perturbation theory provides an 
appropriate description in the weak coupling regime. 




FIG. 1: (color online) Illustration of the correlation energy 
as a function of the coupling strength for the case of N=8 par- 
ticles and Q, = 8. The exact result solution (red solid curve), 
BCS (black dashed line), standard perturbation theory (green 
open squares), Eq. ((3J and QP 2 T theory with second order 
correction Eq. (|14[) (blue filled circles), are displayed. 



In figure[TJ an example of standard perturbation theory 
(SPT) is presented for the TV = 8 particles and constant 
coupling case, i.e. Vij = —g. The correlation energy 
defined as the difference between the Hartree-Fock en- 
ergy Eq and the ground state energy Eo obtained with 
STP are compared to the exact solution and BCS result. 
In the latter case, pairing correlation is non-zero only 
above the threshold value g/Ae ~ 0.3. As illustrated 
from Fig. [TJ standard perturbation theory matches with 
the exact result below the threshold but significantly un- 
derestimates the correlation for larger g value. This as- 
pect underlines the highly non-perturbative nature of the 
pairing quantum phase-transition. On the opposite, one 
of the advantage a theory like BCS is the possibility to 
incorporate non-perturbative physics even in the strong 
interaction case by breaking the U(l) symmetry associ- 
ated to particle number conservation. 



III. QUASI-PARTICLE PERTURBATION 
THEORY 



To provide a proper description of both the weak and 
strong pairing strength regime, it seems quite natural to 
try to combine theories based on quasi-particlcs and pcr- 
turbative approaches. This possibility has been explored 
long time ago in Refs. pa . I2H ] mainly to discuss the 
removal of " dangerous diagram" occuring in normal per- 
turbation theory. Note that recently, it has also been re- 
visited as a p ossible tool to perform ab-initio calculation 
in nuclei [28| based on Gorkov- Green function formalism 
[13 . In the following, it is assumed that the BCS/HFB 
approach has been applied in a preliminary study and 
that the hamiltonian (JT|) is written in the canonical ba- 
sis of the quasi-particle ground state. Then, the ground 
state takes the form 



i*o>= = n(^+^ a K)i->' 



(5) 



i>0 



and is the vacuum of the quasi-particlc creation operators 
defined through: 



0} 



Uia] - Vidz, 
Uiot + V l a l . 



(6) 
(7) 



In the HFB/BCS theory, the original hamiltonian is re- 
placed by an effective Hamiltonian that is conveniently 
written as : 



Hq — Eq 



E, 



(8) 



where Eq is the BCS/HFB ground state energy, while Ei 
corresponds to the quasi-particle energy given by: 



V(e t - A) 2 + 



(9) 



where A is the Lagrange multiplier used to impose the 
average particle number while A is the pairing gap (for 
a detailed discussion see [lj|). 

For even systems, excited states of Hq are 2 quasi- 
particle (2QP), 4 quasi-particle (4QP), ... excitations 
with respect to the ground state. The original hamilto- 
nian H contains many terms that are neg lected in H Q 
and that are responsible from the deviation between the 
quasi-particle and the exact solution. However, noting 
that: 

ff|*o> = (h -^VijUfVftiflp^ |$ ), (10) 

it can be anticipated that the main source of discrepancy 
is due to the coupling of $0 with the 4QP states. Some 



1 Note that here, it is implicitly assumed that H is replaced by 
H — XN 
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arguments showing that 4QP states should improved the 
description of pairing especially in the weak coupling 
regime have been given in ref. [29j . 

Below, perturbation theory is applied assuming that 
Hq (Eq. [8]) is the unperturbed hamiltonian while the 
perturbation V is given by: 



(11) 



V couples the ground state with the 4QP states, defined 
as (for i > j): 



I* 



/3j/3^]/?j|$ ) 



and associated to the unperturbed energy 
Ea =E + 2(E i +E j ), 



(12) 



(13) 



The present approach, that is a direct extension of 
standard perturbation theory, is called hereafter quasi- 
particle perturbation theory (QP 2 T). Using Eq. ([3]), the 
second order correction to the ground state energy is 
equal to: 



E 



(2) 



i>j 



Et + E, 



(14) 



This correction properly extends Eq. (j4]) from the normal 
to the superfluid phase. Indeed at the threshold value of 
g, i.e. when A — > 0, the 4QP states identify with 2p-2h 
excitations while: 



Ei + Ej -> \ei - A| + \sj - A| = Ei - Ej 



(15) 



and the standard perturbation theory case is recovered. 

The result obtained with the QP 2 T approach at sec- 
ond order in perturbation (|14|) arc displayed in figure 
[1] with filled circles. Note that below the BCS thresh- 
old, standard perturbation theory is used. The present 
approach can be regarded as a rather academic exercise 
but it turns out to provide a very simple way to extend 
mean-field theory based on quasi-particlc states. In par- 
ticular, it avoids the threshold problem of the latter and 
improves the description of correlation in the intermedi- 
ate and strong coupling case. 



IV. EFFECT OF THE RESTORATION OF 
PARTICLE NUMBER 

Similarly to the original quasi-particlc theory, the en- 
ergy deduced from the QP 2 T contains spurious contribu- 
tion coming from the fact that the perturbed state does 
not preserves particle number. Indeed, using standard 
formulas, to second order in perturbation, the ground 
state expresses as: 



l*o> 



l^ 0) ) + in 



I* 



(2)\ 



i>j 



(16) 



where I'J'o }, i — 0, 1, ... denotes the contribution to the 
state at i th order in perturbation and where 



{ufv[+u]m 

Ej + Ej 



(17) 



Above the BCS threshold, neither |$o) nor \$i.j) are 
cigenstates of the particle number operator. 

The most direct way to remove spurious contributions 
due to the mixing of different particle number is to in- 
troduce the operator P N 



1 



P* = ^ dip e ^ N - N \ 
that projects onto particle number N. 



(18) 
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FIG. 2: (color online) Illustration of the correlation energy 
as a function of the coupling strength for the case of N=2 
particles and Q — 2. The exact result solution (blue solid 
line), PAV result (green short-dashed line) and the combina- 
tion of projection and perturbation theories (red filled circles) 
are shown. The original BCS (black long dashed line) is also 
shown as a reference. 
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FIG. 3: Same as figure [2] for the N = SI = 8 case. 

The most straightforward way to combine projection 
with quasi-particle perturbation theory is to directly take 
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FIG. 4: (color online) Same as figure [5] for the N = O = 12 
case. 



the expectation value 



En 



(19) 



with l^^} = P N \^o) and where \^o) is truncated at a 
given order in perturbation. This approach will be re- 
ferred as the projected quasi-particle perturbation the- 
ory (QP 3 T) in the following. When only the zero order 
in perturbation is retained in Eq. (|16p . the QP 3 T iden- 
tifies with the projection after variation (PAV) that is 
commonly used, especially in the nuclear Energy Den- 
sity Functional approach (EDF) [2(| ■ Formulas useful 
to compute the expectation values of one- and two-body 
operators with projection are given in appendix [A] The 
projection is performed numerically using these expres- 
sions and the Fomenko discretization procedure of the 
gauge-space integrals [HI, [HI- Here, 199 discretization 
points have been used. Note that the number of points 
can be reduced down to 5 without changing the result. 
The correlation energies obtained for N = 2, N = 8 and 
N = 12 (each time with = N) using the second order 
QP 3 T arc shown in figurcs[2J[3]and[I]respectively. In each 
case, the original BCS, the PAV and the exact solution 
are also shown. The result of QP 3 T almost superimposed 
with the exact solution. Only a slight difference can be 
seen around the BCS threshold. 

Besides the energy, other observables can also be esti- 
mated. As an example, the quantity n,(l— rij) is shown in 
figure[5]as a function of single-particle energies, where rii 
are the occupation numbers of single-particle states. This 
quantity is a measure of the deviation from the indepen- 
dent particle picture where it is strictly zero for all states. 
The use of perturbation and projection considerably im- 
proves the description of one-body observables especially 
at intermediate coupling where the original BCS deviates 
significantly from the exact solution. In addition, even if 
the PAV differs significantly from the exact case as it was 
noted in ref. [IH, the extra mixing with the 4QP states 
compensates this drawback of PAV. 
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5: (color online) Example of evolution of the quantity 
— m) obtained with QP 3 T (filled blue circles) as a func- 
tion of single-particle energy for the N = £7 = 8 case with 
g/Ae = 0.4 (a) and 0.8 (b). In both case, the exact result 
(red solid line) the BCS result (black dashed line) and the 
PAV result (green open squares) are shown. 



In view of this agreement, it seems that the QP 3 T 
does automatically select important many-body states, 
namely projected ground state and projected 4QP states, 
on which the true eigenstatc decomposes. These states 
are highly non-trivial multi-particle multi-hole mixing 
that can also be describe by direct diagonalization of 
the Hamiltonian but that is much more demanding nu- 
merically. Indeed, the size of the 4QP Hilbert space is 
Q(Q — 1)/2 while the size of the matrix to diagonalize the 
hamiltonian is (fi!/[iVp a i r !(f) — V pa j r )!]). It is important 
to recall that here no diagonalization is required since 
the mixing coefficients are directly given by the quasi- 
particle perturbation theory, Eq. ([T7|) . These features 
make the approach rather simple to implement on exist- 
ing HFB/BCS codes to provide a much better approx- 
imation than the PAV that is often currently used. In 
addition, by contrast to the Variation After Projection 
that is rather involved [HI . fl5j , no extra minimization is 
required. 

Contrary to the exact diagonalization, the QP 2 T and 
QP 3 T can be performed even for large particle number. 
As an illustration, in figure [51 a systematic study of the 
correlation energy evolution obtained with some of the 
approaches presented above as the number of particle in- 
creases up to N = 100 for the case g/Ae = 0.8. With 
standard diagonalization techniques, the exact solution 
can hardly be obtained for N > 14. Other approaches 
based on quasi-particle theories can be applied without 



5 



difficulties. Note however, that the QP 3 T requires to per- 
form more and more gauge angles integrations (see ap- 
pendix |A| as N increases making the calculations more 
time consuming with respect to the non-projected theo- 
ries like BCS or QP 2 T. 

From this figure, we also note that differences between 
theories are seen only for rather small particle number 
N < 30, while above BCS and QP 3 T cannot be distin- 
guished. However, as it has been discussed previously, 
for N < 30, the QP 3 T is the only theory that can pro- 
vide an excellent reproduction of the exact result when 
available. 
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FIG. 6: (color online) Evolution of the correlation energy per 
particle as a function of the number of particles for g / Ae = 0.8 
for the BCS (black dashed line), QP 2 T ( green open squares) 
and QP 3 T (blue filled circles). The exact case (solid red line) 
is also shown up to N = 14 particles. 



V. SUMMARY 

In this work an extension of the standard many-body 
theory to treat the pairing problem is introduced. In- 
cluding the effect of 4QP pertubatively to extend usual 
BCS/HFB removes the problem of sharp transition from 
normal to supcrfluid phase and significantly improves the 
description of pairing. Last, when restoration of particle 
number is performed within the perturbative approach, 
a perfect agreement with the exact results is found. This 
finding provide a direct proof of the importance of 4QP 
state to extend mean-field theories. In addition, it is 
shown that the quasi-particle perturbation theory can 
be implemented even for large particle number without 
special difficulties, making the technique rather attrac- 
tive and much simpler than other approaches, like full 
diagonalization, variation after projection or Quantum 
Monte-Carlo techniques. Last, it is worth mentioning 
that the present technique can be directly and rather 
easily implemented on existing BCS/HFB codes (23l. |24| 
to improve the description of pairing correlations. 



Recently, the use of Gorkov-Green function theory has 
been proposed [28| as a possible tool to perform ab-initio 
calculation for nuclei. This theory provides a general 
formalism based on quasi-particle states. The result ob- 
tained in the present study are rather encouraging to 
pursue in that direction and that projection might be 
needed. 
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Appendix A: Expression of projected quantities 

Starting from the standard expression of the quasi- 
particle ground state (Eq. ([5])), the 4 QP states are given 
by: 

i*«> = ( ^ • ' .">:) ( v s ■ i ^y,) 

J] (U k + V k a\a^\-). (Al) 

For compactness, this expression is written as 

l*m> = H(u? + Vral4)\-). (A2) 

k>0 

This notation includes the ground state case (to = 0). 
The state obtained in QP 3 T can be generically written 
as \^q) = Y^m Cro-Pjvl^m)) and for any operator O that 
conserves the particle number, we have 

=J2 C n C ™(*n\OP N \$ m ) (A3) 

where 

/■27T p—iipN 

{® n \OP N \$ m )= J dtp-^-(* n \0\$ m (<p)), (A4) 
and 

i*i»(v)> = n(^ m+y ™ e2iva H)i->- ( A5 ) 

fe>0 
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Starting from this expression it could be deduced that 

"2tt 



($ n \P N \$m) = y I II ( U k U ™ + V k n V k m e 2i *) (A6) 

77 J ° k>0 



(#„ \a\ ai P N |$ m ) = i- d V e-^ N VrV % m e 21 * J] (^C/f + V^e 2 * 4 



Z7T JO ,s 



fc>0,fc^(i,j) 

where the latter expression is valid for i =/= j. 
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